杨昕衍:
library(ChIPpeakAnno)
杨昕衍:
OCT4.DP <- toGRanges("OnlyDP/DP.Venn_OCT4motif.filter.bed", format="BED", header=FALSE)
SOX2.DP <- toGRanges("OnlyDP/DP.Venn_SOX2motif.filter.bed", format="BED", header=FALSE)
SOX17.DP <- toGRanges("OnlyDP/DP.Venn_SOX17motif.filter.bed", format="BED", header=FALSE)
SOX15.DP <- toGRanges("OnlyDP/DP.Venn_SOX15motif.filter.bed", format="BED", header=FALSE)
AP2.DP <- toGRanges("OnlyDP/DP.Venn_AP2motif.filter.bed", format="BED", header=FALSE)

ol.DP <- findOverlapsOfPeaks(OCT4.DP, SOX2.DP, SOX15.DP, SOX17.DP, AP2.DP)
杨昕衍:
makeVennDiagram(ol.DP, fill=c("#8C9DD5", "#E68D8D", "#8CCF8C", "#F5F4A0"))
杨昕衍:
# I need to extract 48\33\224\36\19
S15_O4_AP2.DP.Venn <- as.data.frame(ol.DP$peaklist$`OCT4.DP///SOX15.DP///AP2.DP`)
S15_O4_AP2.DP.Venn <- S15_O4_AP2.DP.Venn[,1:3]
write.table(S15_O4_AP2.DP.Venn, "S15_O4_AP2.DP.Venn.bed", row.names = F, col.names = F, quote = F, sep = "\t")



bedTools.2in<-function(functionstring="bedIntersect",bed1,bed2,opt.string="")
{
  #create temp files
  a.file=tempfile()
  b.file=tempfile()
  out   =tempfile()
  options(scipen =99) # not to use scientific notation when writing out
 
  #write bed formatted dataframes to tempfile
  write.table(bed1,file=a.file,quote=F,sep="\t",col.names=F,row.names=F)
  write.table(bed2,file=b.file,quote=F,sep="\t",col.names=F,row.names=F)
 
  # create the command string and call the command using system()
  command=paste(functionstring,"-a",a.file,"-b",b.file,opt.string,">",out,sep=" ")
  cat(command,"\n")
  try(system(command))
 
  res=read.table(out,header=F)
  unlink(a.file);unlink(b.file);unlink(out)
  return(res)
}

bedTools.2in("bedIntersect",bed1,bed2) # equivalent to "bedIntersect -a bed1 -b bed2" on terminal


intersectBedFiles.foverlaps <- function(bed1,bed2) {
  require(data.table)
  bedKey <- c("CHROM","START","STOP")
  if(nrow(bed1)>nrow(bed2)) {
    bed <- foverlaps(bed1, bed2, nomatch = 0)
  } else {
    bed <- foverlaps(bed2, bed1, nomatch = 0)
  }
  bed[, START := pmax(START, i.START)]
  bed[, STOP := pmin(STOP, i.STOP)]
  bed[, `:=`(i.START = NULL, i.STOP = NULL)]
  if(!identical(key(bed),bedKey)) setkeyv(bed,bedKey)
  if(any(bed[, STOP+1 >= rowShift(START), by=CHROM][,V1], na.rm = T)) {
    bed <- reduceBed.GenomicRanges(bed)
  }
  return(bed)
}

rowShift <- function(x, shiftLen = 1L) {
  #Note this function was described in this thread:
  #http://stackoverflow.com/questions/14689424/use-a-value-from-the-previous-row-in-an-r-data-table-calculation
  r <- (1L + shiftLen):(length(x) + shiftLen)
  r[r<1] <- NA
  return(x[r])
}

reduceBed.GenomicRanges <- function(bed) {
  setnames(bed,colnames(bed),bedKey)
  if(!identical(key(bed),bedKey)) setkeyv(bed,bedKey)
  grBed <- makeGRangesFromDataFrame(bed,
    seqnames.field = "CHROM",start.field="START",end.field="STOP")
  grBed <- reduce(grBed)
  grBed <- data.table(
    CHROM=as.character(seqnames(grBed)),
    START=start(grBed),
    STOP=end(grBed),
    key = c("CHROM","START","STOP"))
  return(grBed)
}